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ABSTRACT 

In this paper we present theoretical and numerical results for inverse 
problems involving estimation of spatially varying parameters such as 
stiffness and damping in distributed models for elastic structures such as 
Euler-Bernoulli beams. An outline of algorithms we have used and a summary of 
our computational experiences are presented. 
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1. Introduction 


In the past several years there has been an increased interest in 
the use of continuum models [13], [18], [22] , [27] , [28], [29], [30], 
in the study of large complex connected structures such as those being 
planned for deployment in space. These flexible platform and antenna 
structures, which are frequently composed of large lattice, panel and/or 
beam-like components, are typically constructed of graphite epoxy 
composite materials. Preliminary experimental testing suggests that one 
can expect significant material property changes (e.g., decreases in 
material damping by as much as 500%) due to ageing, environmental stress, 
fatigue, etc., during periods of deployment of structures composed of 
these composite materials. Therefore the identification or estimation of 
structural parameters (e.g. , bending and shear rigidity, moments, damping, 
loading) will play an important role in the modeling , control and 
stabilization of these large space structures. 

In recent efforts we (along with some of our colleagues and associates) 
have contributed to the research literature [2] , [3] , [4] , [5] , [8] , [10] , 
[11] , [19] , [25] in developing methodology (theoretical as well as computa- 
tional) for such problems. In this paper we detail and extend some 
earlier results reported in [2], [4], [5] and [10]. Our specific aims here 
are twofold: (i) to present some of the ideas behind the theoretical 

results stated in [4] and [10] , and (ii) to present a more extensive 
summary of some of our numerical findings for estimation of spatially 
dependent parameters, especially damping (again, earlier findings were 
presented in [4] , [5] and [10]) . To be more precise, we note that in 
many of the earlier efforts cited above, the focus was on convergence 
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results in a functional analytic framework (e.g. , semigroups, dissipative 
operators, sesquilinear forms) and entailed certain smoothness (on the 
approximation elements) and compactness (on the admissible parameter 
sets) hypotheses. Here we present in §2 a theoretical approach (presented 
earlier in [4] , [5] and treated in the context of hybrid systems describ- 
ing the undamped vibration of beams with tip appendages in [10]) based 
on weak or variational arguments in the spirit of those from 
widely known "finite-elements" approaches to the approximation of partial 
differential equations (for a summary and numerous literature references, 
see [14]). This approach permits, in addition to relatively weak smooth- 
ness assumptions on approximation elements, a substantial relaxation of 
compactness assumptions on the admissible parameter sets. We comment 
further on this aspect of our presentation in §5 below and also refer 
the reader to [1] for a more complete discussion of weak vs. strong 
formulations in the context of inverse or parameter estimation problems. 

In §3 we discuss implementation of some approximation ideas as we 
have employed them. Included is brief mention of the three principal 
optimization schemes we have used (sometimes in a hybrid method) : a 

finite-difference Levenberg-Marquardt algorithm, a conjugate-gradient 
algorithm, and the popular Broyden-Fletcher-Goldfarb-Shanno algorithm. 

A summary of some of our numerical findings for a particular model (a 
damped Euler-Bernoulli beam) is given in §4 while the last section is 
devoted to a brief discussion of related and continuing efforts. 
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2. Theoretical Considerations 

In this section we give a brief outline of convergence arguments 
associated with the algorithms and approximation ideas that we have 
used in obtaining the numerical results presented in subsequent sec- 
tions. The arguments are applicable to quite general models of practical 
interest (e.g. , see [10]) and we choose a specific simple model 
- a cantilevered Euler-Bernoulli beam with viscoelastic damping - only 
to illustrate the theoretical ideas. 

We assume a normalized (mass = 1, length = 1) thin elastic beam with 
built-in or clamped end at x = 0, free end at x = 1, which is subject 
to Kelvin-Voigt damping. The equations for transverse vibrations 
embodied in the Euler-Bernoulli theory are given by (here D = 9/9x) 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 

(2.4) 


U tt + ° 2 { <!i d2u + vV = f ' 

u(t,0) = Du ( t , 0) = 0, 
{q 1 D 2 u(t,*) + q 2 D 2 u t (t,*)} x _ 
[D{q 1 D 2 u(t, • ) + q 2 D 2 u t (t,‘) }] 


0 < x < 


t > 0, 



1 , 


with initial data 


(2.5) u (0) = u Q 

(2.6) u t (0) = v Q . 

We assume here, for ease in exposition, that the initial data u^, 
v Q are independent of the parameters q = (q^q^. The case u Q = 
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Ug(q) / = v Q (q) can be treated with ideas similar to those given 

below. In the general case the loading term f may also depend on 
parameters to be estimated; again the arguments below are easily 
extended to cover such situations and we consider only the simple 
case f = f (t ,x) . 

The parameters q are restricted to lie in the admissible set 
Q = (q = (q 1 /q 2 ) J q ± € c[o,i] , o < 1 q i (x) £ v} 

where C \' C 2 ,V are ^•'- xe ^ constants. We remark that if no damping 
is present (q^ = 0) , the convergence arguments outlined below can 
be made with even less effort. 

We reformulate the system .(2.1) -(2. 6) in weak or variational 
form in the state spaces V and H = L 2 (0,1) where 

V = H* = O € H 2 (0,1)|>H0) = DiMO) = 0} . 

Denoting the usual inner product in L 2 (0,1) by we re- 

place (2.1)— (2.4) by 

(2.7) + <q 1 D 2 u + q 2 D 2 u fc ,D 2 \p^> = <f,i 

2 

and seek solutions u with u(t) £ satisfying (2,5) (2.6) and 

(2.7) for all ip £ H 2 . 

Remark ; As we have noted, other models can also be readily treated 
with the approach given here. For example, in the case of a simply 

supported beam, the boundary conditions (2.2) -(2.4) are replaced by 

2 2 1 
u(t,g) = D u(t,n) = 0, n = 0,1 and the state space V = H (0,1) ft H Q (0,1) 

2 

is used in place of H*. 
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When for example a tip mass of magnitude m is rigidly attached 
to the free end of the beam, the natural boundary condition (2.4) 
expressing zero shear is replaced by the ordinary differential 
equation 


mu tt (tf1 ^ " l E> {‘i 1 D2u ( t '*) + q 2 o 2 u t (t, •) }] 


= g(t) 


X=1 


where g is an external load applied to the tip mass in the transverse 
direction. In this case the appropriate choice for H is the product 
space R x L^CO,!) with 


v = Un,40 € H : ijj 6 h ( 0 , 1 ), HO) = DtMO) = o, n = Ml) } 


(see [10]) . More will be said about our general approach in the 
context of hybrid systems (i.e. coupled systems of ordinary and 
partial differential equations) when we discuss specific examples 
and our numerical findings in Section 4. 

It is not difficult to use standard arguments to show that (2.5), 

(2.6), (2.7) is well-posed. That is, under reasonable smoothness 

assumptions on f, u^, q^, q^ and the positivity constraints 

of Q on q^, q 2 , one can modify the arguments in [20, p. 272-281] 

to obtain existence of a unique solution u (on any finite interval 

[0 ,T] ) satisfying u G C ( [0,T] ,V) , u fc € C ( [0 ,T] ,H) , u G L 2 ([0,T],V') 
2 

where V = H*(0,1) , H = (0,1) and V' is the dual of V with H 

as pivot space. Under additional smoothness hypotheses, one can 
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argue that this solution is actually a strong solution of (2.1)- 
(2.6) enjoying stronger smoothness properties (for example, see the 
results given in [10] , [11]). 

We consider a class of problems for the estimation of the 
parameters q^q^ 9-^ ven observations of the system (2.1) — (2.6) . 
Given observations u^ for u(t.,x.) and a parameter set g eg, 
we seek q* € Q to minimize over Q the least-squares criterion 


( 2 . 8 ) 


J (q) 


l |u(t ± ,x ,;q) 

i/ j 



2 


where u = u(q) is the solution to (2.5)— (2.7) corresponding to q. 

Without additional assumptions on Q, these problems are in- 
tractable (from both a mathematical and computational viewpoint) . We 

2 

shall assume throughout that Q is compact in the = C([0,1] ,R ) 
topology. Even with this compactness assumption, the minimization 
problem for (2.8) is infinite dimensional in both state and parameter 
space and thus is not readily solved without approximations . 


N 2 , N 

For state space approximations H <= H* we let H be finite 

N N 

dimensional, N = 1,2,..., and let P:L(0,l)-*-H be the orthogonal 

N N 

projection of L^(0,1) onto H . We assume that H satisfies: 

2 N 2 

(2.9) For each ip € 1^(0, T;H*) , P 4> ip in (0,T;H^) . 


These hypotheses are satisfied by a number of useful and popular 

families of approximations: quintic or cubic B-splines modified to 

2 

satisfy the boundary conditions defining H* (see [ 2 ] , [26]) 

Hermite cubic splines (modified); etc. 
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M M 

For parameter set approximations Q , we suppose that Q is 
M M 

given by Q = where the set Q and the mapping : 

2 

Q-> C ( [0,1] ,R ) satisfy: 


(2.10) Q M is compact in the So topology; 


(2.11) i^ (q) q in the SZ topology, uniformly in q £ Q. 


These conditions are relatively mild and include some practically 
useful approximation schemes such as linear and cubic interpolatory 
splines (for discussions, see [7 ] , [ 9 1) • We note that we do not 


_M 


require Q c Q, although in some situations this may be automatically 


satisfied. In other cases it may be desirable to impose the constraints 

M 

in Q explicitly in using the sets Q in computational examples. 

N 

For any q 6 Q, we may define in H approximating systems for 

N N 

(2.5)- (2.7) as follows: we seek u (t) £ H such that for all 


\p £ H 


N 


^ 2 ^ 2 2 
(2.12) + < ^ q i D u + q 2° u t ,D = 


(2.13) u N (0) = P^u 


N N 

(2.14) u t (0) = P v Q 


We then define, for approximation indices (N ,M) , the approximating 
estimation problems: 
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Minimize 


(2.15) 


N r i N ~ \ 

J (q) = I |u (t^x.jq) - u | 

i/ j 3 3 


M N 

over q £ Q , where u is the solution of (2. 12) - (2 . 14) . 

-N 

Let q be solutions of the (N,M) estimation problems, N = 1,2, 
M 

..., M = 1,2,... (Such solutions exist since Q M is compact in $£ 

and q J N (q) is continuous in the 5# topology.) Since q^J £ Q M = 

— N — N ~N , ^ 

i (Q), there exist q in Q such that i.,(q , ) = q, . The j& com- 
M M M M M 


pactness of Q implies the existence of a convergent subsequence 
r_N i r i N . 

jq M 3 j of |q M | With q^ -> q £ Q as N j + 08 , \ * 00 • 

The limit function q is an obvious candidate for a solution to 


the problem for (2.8) . To see that it does indeed provide a minimum 

for (2.8), we first observe that property (2.11) for i along 

M 

with the inequality 


q 3 - q 

M, M 

< 

k 



N. 




N. 


N. 

) - q 3 

+ 

— T 

V " q 

k 


k 


r N o 

k) 


guarantees that 
definition 

N. / N.v N. 

j 3 k) - j 3 


\ -k 

A 

also converges to q in jf- Next we have by 


M. 


(q) for all q £ Q 


or 


N. / N.v N. 

J 1 J ki M (q>> for all q £ Q. 


V 


M, 


Thus, taking the limit as °° yields the desired in- 

equality 
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J (q) <_ J (q) for all q € Q, 


N n 

if we can argue J (q ) -*■ J (q) as N -* 00 , n -»■ 00 , for any sequence 

{q 11 } with q 1 -*■ q in But this follows immediately once we have 

argued that u N (t.,x.;q n ) -*• u(t.,x.;q) for arbitrary q 11 q. The 

. J i 1 

remainder of this section will be devoted to a sketch of arguments 
for this convergence. 

We recall (2.5)-(2.7) and (2.12) - (2 . 14) where we must consider 

N, n, , . n . . 

u (q ) , u(q) with N *■> °°, n -* °°, q -> q. Simple remdexing argu- 

N N 

ments reveal that.it suffices to argue that u (q ) -*■ u(q) as 

N N 

N -*■ “ for arbitrary sequences {q } with q -*■ q as N -> “. Here, 
N 

of course, u (q) and u(q) are the solutions of (2 . 12) - (2. 14) and 

(2.5)— (2.7) , respectively, corresponding to q. 

N N N N 

Let q -*• q be arbitrary and let u ,u denote u (q ) ,u(q) 

throughout below. We define 


N N 
z = u 


, N x N . 

(q ) - P u(q) . 


Then (2.5), (2.6), (2.13), (2.14) imply that z N (0) = z^(0) = 0 (in 

cases where u q ,v q depend on q, the arguments differ slightly since 
N N 

then z (0) ,z^(0) are not zero, but approach zero as N “ under 

suitable smoothness assumptions on u Q ,v 0 ) . Using (2.7) and (2.12) 

N 2 

we have for G H c: h a 
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/ N . V ^ N , N . V 

<Z tt >> = <U tt “ U tt + U tt " P U tt A > 

. 2 N 2 N 2 \ ^ 2 N 2 N 2 v. 

(2.16) = <q.jD u - q x D u , D \p/> + <s,q 2 D u fc - q 2 D u fc ,D \p 

+ <(I - P N ) U . /^> - 

TIT- 

(If f depends on unknown parameters q we would also have a term 
<^f (q N ) -f (q) above but the essential features of our presentation 

would again remain the same.) 

Adding appropriate terms to both sides of (2.16), we find 


^ N , \ / N 2 N 2 , \ / n z N 

< C z tt '4'> + z ,d i)j> + <q 2 o z t ,D ijj> 

/ 2 N 2 N 2 n. ,< 2 N2N 2 ^ 

= <^0 u - q^D P u,D <J>> + <q 2 D u fc - q 2 D P u t> D 

+ <(I - P N )u tt ,^> . 


N 2 N 2, 


(2.17) 


N 


N. 


Choosing i/» = z (which is in H ) , we obtain 




✓ .N 2 N\ / ,N Z N ** , IN N \ 

= <A 1 ,D z t > + < A 2t /D z fc > + <5 , z t > 


N _2 N • 


N _2 N 
)l q 2 ° Z t 


,N N 


N _ 2 N 2 N , ,N _ N 

where A. =q.Du-q.DPu and o = (I - P )u , . 

x x x tt 


This implies 

(since 

q 2 (x) c 2 ) 




2 



I Aj ! 

N 

_i_ 

N „2 N 1 

J% D 2 1 r °2 

r, 2 N 

D z t 

2 dt l 

zl* 

t 

T 

■ 


d ^ . N 2 N — ^ LN Z IM\ 

±dt < A l'° 2 < A lt ' D z > 

✓ „ N N\ 
+ <6 ,z fc > 


N 2 N« 


4c, 



2 

+ c 

2 N 
D z 

2t 

2 

t 


2 
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(If damping is present, but we are only guaranteed q 2 (x) — 

the arguments must be modified at this step, but again the essential 

ideas are similar to our presentation here.) From this last inequality 
we have 


d 

dt 



+ 



2 N 
D z 



2 N 
z 



< 2 





2 


+ 2<S N 


N 

z 

t 


>. 


Integration of this inequality leads to 


/ N 2 N 

k/ q i D 2 


N , 

Z t (0) 


IP" 


2 N,„, 
z (0) 


+ 2<A^,D 2 z N > 


<A X (0) , D Z (0)> 


ft { ; 

n *• 


/.N N\ 

<A lt ,D z >| 


2c, 


2t 


+ 2 




In the case we are considering here (u Q ,v 0 independent of q) , we have 

N 2 N 

z^{0) = 0 and one also can easily argue that D z (0) = 0. We there- 
fore find 


N 

2 

2 N 

2 < j_ 

A N 

2 c 

1 

2 N 

z 

+ c. 

D z 


A_ 

+ ~ 

D z 

t 

1 


" c i 

1 

2 




ds. 


Defining 
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v N (t) 



and 



we have obtained 

N N f fc N 

v (t) F (t) + v (s)ds . 

■'O 

N 

An application of the Gronwall inequality yields v (t) -* 0 if we 

N N 2 

can argue F (t) -*■ 0 for each t. But since z (t) 6 H* , the con- 

N N 2 . 

vergence v (t) 0 actually yields z (t) -* 0 in H* which in 

turn yields the desired convergence. 

N 

To complete the arguments, one recalls the definition of 6 

and A^. Under assumptions (2.9) and the compactness of Q in ^ 

(i.e., q N •* q in , the arguments for F N (t) -*■ 0 can readily be 

reduced to smoothness requirements on u = u(q) (e.g., u,u fc £ 

2 

L (0 ,T; H ) , u £ L_(0,T; L-.) ) - These in turn can be established as 
2 * tt 2 2 

indicated above under additional smoothness hypotheses on the data 
(e.g., UQ,v Q ,q,f) in our problem. 
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3 . Implementation and Computational Considerations . 

We next discuss implementation of schemes for the sequence of 
approximate estimation problems involving the systems (2.12) - (2.14) 
and criteria (2.15). This involves Ritz-Galerkin type procedures 

embedded in optimization routines. For these one needs state approxi- 

N M 

mation subspaces H and parameter approximation sets Q . In the 

N 

case of the work reported on m this paper, we employed subspaces H 

generated by either cubic or quintic B-splines appropriately modified 

to satisfy the essential or geometric boundary conditions (2.2) of 
2 

V = H*(0,1) (or the analogs of these conditions in the event we are 

considering some other type of boundary condition, e.g. simply- 

supported, or tip appendage) . We have used either linear or cubic 

M 

interpolatory splines to define the sets Q . Once approximation sets 

and have been chosen the equations (2.12) reduce to matrix differen- 

N 

tial equations for the "Fourier" coefficients of u relative to the basis 

N M 

elements for H and Q . 

To illustrate the ideas in a specific case, we consider cubic spline 
state approximations and linear spline parameter approximation sets for a 
cantilevered damped beam modeled by (2.1)-(2.6). 

N 

We first describe construction of the basis elements for H . For 

N r lK 

a positive integer N and partition A = {x^K_ Q , = i/N , of [0,1] f 

3 N N 

we denote by S (A ) the set of cubic splines with knots A (i.e., 

the set of functions s such that s is a cubic polynomial on each 

interval and is c on [0,1]). Let lB^.L_ ^ be the 

3 N 

standard cubic B-spline basis set for S (A ) - see, for example 
[24, p. 79]. The cubic B-spline B^ has support in ^ x i_2 ,x i+2^' 
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has values 1, 4, 1 and slopes 1/N, 0, -1/N at the knots x^ x^, x i+1 

3 N . . 2 

respectively. The set S (A ) is not contained in H* but we use the 

-N N N N 

elements B. to generate basis elements B., that satisfy B.(0) = DB.(O) = 
3 3 -33 

0. More precisely, we define 


N -N „ ~N N 
B l = B o - 2B i - 2B -l 


N -'N 

B ± = B ± , i = 2,3,. . . ,N+1 


, , N . N N N , 

and take H = span{B, , B_ , . . . ,B„ , , ; . 

12 N+l 


Remark : When considering the clamped- free beam with tip mass, we choose 

H N = span (3^, 3^ 3^+1 } 


N N N 

3. = (B. (1) , B.) , j = 1,2,... , N+l . 

3 3 3 

M M 

To construct Q for a positive integer M, we let L(A ) be the 

set of piecewise linear splines [26, p.10; 24, p. 48] corresponding 

M . M 

to the partition A = {i/MK . Basis elements for this set are given 

M M 

by the standard "hat" functions b^, i = 0,1,..., M where has 

support in (x^ ^,x^ + ^) with values 0, 1, 0 at x i+l re “ 

spectively. Letting i denote the usual interpolation operator for 
M 

L(A ) - see [26, p. 10] , we then define for a given Q, the approxima- 

M 

tion set Q = i (Q) . Note that in this particular case, the constraints 
M 

M 

c^ q^(x) v are preserved m Q . 

We thus have that any solution u N (t) 6 h N of (2.12) has the 
representation 

N+l 

(3.1) u N (t,x) = £ w^(t)B N (x) 

3=1 D 3 

M 

while the coefficients q^,q 2 to be chosen from Q have the repre- 


sentation 



-15- 


M V M 

(3.2) q"(x) = l q b"(x). 


j=0 


For these fixed bases (for given indices (N,M)), it is easy to 

argue that (2.12) reduces to a second order N+l dimensional system 
~N N N T 

for w = (w^, . . . ,w N+1 ) . Specifically, (2.12) must hold for each 
N 

ip - B , k = 1, . . . ,N+lj this, in view of the representation (3.1), yields 
an N+l system of second order ordinary differential equations which 
can be rewritten as a 2N+2 first order system 


(3.3) 


_N* N, , N N, . , _N 

QW (t) = K W (t) + f 


Q W (0) = W Q 


N N NT N - *N 

where now w = (w /W. 0 ) with w = w., i = 1,2 , — ,N+1, 

1 ZN+2 1+N+l 1 

and 


q n = 


N 




NN 
B.B . 
i 3 


k n = 


K l (q i } 


N 

Vq 2 ) 


[Ki(q,)]i, ; 


2_N_2_N 
% D B.D B . , 


f N = 


,N 




rf, 


[ 4 - 


rl 


_N 


■ 0 

u 0 B j ' 

j = 

,1 


N 



V„B . , , 

j = 

•*0 

0 J-N-l 



j = 1,...,N+1, 
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The approximate estimation problems (for given integers N,M) 
thus reduce to minimization of 


(3.4) 


*N,M ~ 
$ (q) 


“ I 


1.3 


I N /fc M, 

! u (t^Xjjq ) - u 


subject to (.3.1), (3.2), (3.3) where now q is the vector of parameters 
q^_. in (3.2). For each set of fixed indices N,M, this problem can 
be successfully treated with a number of different techniques. We pro- 
ceed to outline some of those which we have used. First, however, we 
note that computation of $ in any optimization routine requires 

solution of the semi-discrete Galerkin equations (3.3). This can be 
accomplished by using the Hindmarsh adaptation of the Gear algorithm 
[17]. For problems where q 2 ( x ) = 0 (no damping), the equations can 
be integrated efficiently using the Adams methods which are part of 
that scheme; however when damping is present the equations are stiff 
and we found it necessary to use the routines for stiff systems that 

are part of the Gear algorithm. 

N M 

The evaluation of $ ' is the expensive part of our algorithms, 
involving the integration of moderately stiff systems of 2N+2 dimen- 
sional differential equations. Some computational savings can be achieved 

due to the special structure of the matrices in (3.3): sparse, banded, 

N N 

with symmetry in the matrices and K^. A Cholesky algorithm can 

be used at each integration step to solve for the N+l subsystem involving 
N N N 

Q 2 . Elements of Q and K can be readily evaluated usxng a composite 
two-point Gaussian quadrature. 
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We have used, in the. results reported herein and in earlier efforts 

[5] , three different methods (sometimes in a hybrid scheme involving two 

. N,M 

of them) to carry out the minimization for $ in (3.4) : (i) a finite- 

difference Levenberg-Marquardt (FDLM) algorithm; (ii) a conjugate-gradient 

(CG) algorithm; and (iii) the scheme due to Broyden-Fletcher-Goldfarb- 

Shanno (BFGS) . We describe each briefly below, referring the reader to 

the references given below for further details. We first note that (i) 

N M 

requires only evaluation of $ since derivative information is obtained 

by finite differences; methods (ii) and (iii) require gradients which. we 

« 

have computed using a costate equation approach (explained below) . 


(i) Finite-difference Levenberg-Marquardt : This quasi-Gauss-Newton [23] 

method is designed especially for minimization of least-squares criteria. 

K 

Rewriting $ N,M (q) of (3.4) as l e ? (q) 2 where e(q) = (e, (q) , . . . ,e (q) ) 

1=1 1 K 

is the vector of residuals or pointwise errors, the L-M scheme generates a 
sequence of iterates 



with the scalar y^ chosen to insure positive definiteness of the system 

~ (k) 

matrix in (3.6) and so that d is a descent direction. Here the matrix 
. J k is the Jacobian matrix of e(q) evaluated at q . For y = o 

.K 

the directions thus obtained agree with those of the Gauss-Newton method 
while for y^ large the directions approach those of steepest descent. 
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This method generates descent directions and exhibits rapid 

convergence (superlinear) . However, to approximate the Jacobian in 

(3.6) by finite differences (which is the standard approach for this 

method), the differential equations (3.3) must be integrated p times, 

vhere p is the number of parameters in the vector q = (q ) Thus 

ij 

at each iteration of (3.5), the approximating differential equations 
must be solved p times to generate one descent direction. This is 
not a significant drawback if p is not too large. For spatially 
varying stiffness and damping coefficients q , q 2 in (2.1), the total 
number of parameters needed to obtain reasonable approximations can 
become large. In this event alternate methods may be superior. 

One can avoid methods based on finite-difference gradients if one 
is willing to compute the necessary gradients (e.g. , by variational 
equations or costate methods) and supply them directly to the iterative 
algorithm. While there is a version of the Levenberg-Marquardt 
algorithm (in MINPACK's LMSTR1) which allows user calculated gradients, 
we have chosen to investigate use of other popular algorithms in some of 
the calculations we have carried out. Since the use of linear variational 
equations in calculating the gradients involves solving systems that grow 
in dimension with the number of sought-after parameters (p in the discus- 
sion above) , we have chosen to use a costate formulation in computing 
gradients. In this case only two differential equations (with dimension 
independent of the number of parameters to be estimated) must be inte- 
grated per iteration. The associated iterative algorithms we have used 
are the conjugate-gradient and the BFGS. 
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(ii) Conjugate-Gradient : The CG algorithm is quite well-known (see 

[16, p. 133-136; 15 , p. 91-98] for more details) and we use the for- 
mulation proposed by Nazareth in [21] . Once again the iterates are 
given by (3.5) but the directions are generated by the recurrence 
relationships 



d Q) = -g (1) 

d (k+1) = _ y (k) + (y (k-l)T y (k) /y (k-l)T d (k-l) )d (k-l) 

+ (y (k)T y (k) /y (k ) T d (k) )d (k) 

. (k) (k+1) (k) ... (k) fc . .. ^ , _ ,N,M 

where y = g - g with g the gradient vector (of $ 

- - (k) 

with respect to the parameters q) evaluated at q 

While CG methods exhibit finite-step convergence when the cost 
criterion is a quadratic functional , in general problems such as those 
under investigation here convergence is often slow, particularly in the 
neighborhood of the extremal point. Newton methods are more suitable 
near the extremal and hence it is often advantageous to switch to a 
quasi-Newton method such as the BFGS. 


(iii) Broyden-Fletcher-Goldfarb-Shanno : In this method one again uses 

~ (k) 

the iterative formula (3.5) with the search directions d now given 
by 


5 <k > = -H k g (k > 


(k) 

where g again denotes the gradient and 
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H k+1 = H k + (1 +yVy /(S^Y )5 6^/6^Y 
k k k k k k k k 


- (V^hV^/S^ ' 


* - ~ (k+1) - (k) _ (k+1) (k) 

6 k = q -q ' Y k~ g - q 


H° = I. 


This quasi-Newton method employs "BFGS updates" H for a matrix which 
is an approximation to the inverse Hessian matrix. In practice, the 
method exhibits rapid convergence for the type of problems we have con- 
sidered, particularly when compared to the CG method (ii) in the neigh- 
borhood of an extremum. Second order information is used (i.e., an ap- 
proximate inverse Hessian matrix) even though only gradient computations 
are required. In our implementations, the gradient is computed using a 
costate formulation and hence as in the case of the CG method, only two 
differential equation solutions (state and costate) are required per step, 
as compared to the FDLM (or any other finite-difference based gradient 
algorithm) which requires p (= dimension of unknown parameter vector) 
differential equations be solved per step. 

The step parameter A is determined by a one-dimensional line 

K 

search, with the BFGS and CG methods offering the advantage that exact 
line searches are not necessary. For a further discussion of the BFGS 
method, the reader may consult [15, p. 38-60], 

The method of using costates to compute the gradient in optimization 
problems with differential equation constraints is well-known to investi- 
gators in control and variational theory. These ideas in a form appropriate 
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for parameter estimation problems are described in some detail by 
Chavent [12] for the case of continuous time observation problems. 
A similar formulation for discrete time criteria such as (3.4) can 
be given. We refer the reader to the Appendix where it is shown 
that the desired gradients are given by 


3* N ' M 

3q. . 
ID 


(q) = - p T (t)J^ i ,w N (t)dt 


with 


^ lj J 



where is the N x N matrix with elements 



rl 

'0 


, M 2 N 2 N 
b . D B . D B. 
3 i k 


dx 


and p satisfies an appropriately defined costate (to (3.3)) equation 
(see (A. 7) in the Appendix). 
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4. Numerical Results and Examples 

We have developed and tested software packages based upon both 

quintic and cubic spline state approximations and cubic and linear spline 

parameter approximations as noted above (some of 1 our findings were 

reported in an earlier version of this paper [5]). Numerical testing was 

carried out on a CDC Cyber 173 at NASA Langley Research Center, a Burroughs 

6900 at the USAF Academy, and an IBM 3081 at the University of Southern 

California. The IMSL version of the Gear/Hindmarsh algorithm (DGEAR) was 

used to integrate the approximating differential equations (3.3). One of 

several methods was used to solve the approximating optimization problems: 

the IMSL implementation ZXSSQ of the Levenberg-Marquardt algorithm, the 

BFGS scheme or a hybrid CG/BFGS scheme (i.e., initially using the CG 

algorithm and switching to the BFGS as we neared a minimizer) . 

To test the approximation ideas, the inverse procedure, and software, 

synthetic data was generated for a number of examples in the following 

* * 


manner: "true" parameter functions q , q 2 were chosen and an algorithm 

employing either finite differences or a high order spline based Galerkin 
scheme was used to generate the corresponding solution values u_ = 
u(t^,x_.;q ) at certain points (t^,x^) on a grid. These values (which 
obviously already contain some "noise") were subsequently used as observa- 
tions or input for the estimation algorithms and software being tested. 

A number of test examples with different shaped stiffness and damping 
functions were investigated in this manner. We report in this section 
some of these results. The results we summarize are typical of our 
numerical findings. Each of the particular examples detailed below was 
run using the ZXSSQ package, although we have also enjoyed success with 



-23- 


the other optimization methods discussed above (e.g. , see [5] where some 

of the results given were obtained using the CG and BFGS schemes) . 

* ★ 

In Examples 1 and 2 below, for "true" parameter functions q^, q^ , 

"data" was generated using a finite difference scheme for observation 
points (t ,xj with t^ = i/10, x^ = j/10, i,j = 1,2,... ,10. The test 
model equation used is given by (2.1) - i.e., 

2 2 2 

(4.1) u + D {q^D u + q 2 D u fc } = f, 0 < x < 1, 

with initial data u(0,x) = u t (0,x) = 0 and uniform loading f(t,x) = 10. 

The boundary conditions used were for a simply supported beam, that is, 

(4.2) u(t,0) = u(t,l) = D^u(t ,0) = D^u(t,l) = 0. 

In the results reported below, the index N will always refer to the 
state approximation level with either N+l quintic basis elements (note 
this is a strong formulation with all four conditions of (4.2) imposed 
on the basis elements) or N+l cubic basis elements (a weak formulation 
in which the zero moment conditions are treated as natural boundary 
conditions) . Similarly the index M will denote the parameter approxima- 
tion index with M+l linear splines or M+3 cubic splines in a basis set. 

In carrying out the optimization step, it was often necessary to impose 
the pointwise positivity constraint q^x) rl c 1 > 0 of Q on the stiffness 
estimate q^' generated by the various schemes. We denote by q N,M the 
converged parameter estimates (to be compared with q?, the true parameter 
functions) and tabulate one or more of the, following measures of performance: 
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(4.3) 

$ N ' M 3 

l 1 

i/ j 

N a N,M. ~ | 2 

u (t if x.;q ) - u^\ 

(4.4) 

e N ' M e 

X 

kj - 

-N,Mi 
q i 1 

where 

H in 

(4.4) 

denotes the norm in L (0,1) 


* * 

Example 1 . We consider an undamped beam (q 2 E 0) with stiffness q^Cx) = 

.15 + .10 tanh [5 (x- . 5) ] . We estimate q^ from start-up value q^ = .15 

using quintic spline state approximations, and cubic spline parameter 

approximations. In Figure 4.1 we present a graphical record of the 
~N M 

convergence q^ -> q* while in Table 4.1 we list the corresponding values 
-m M N,M 

for $ and E . In all of our graphs, the true parameter functions 

~N,M 

are represented with dashed cuives while the estimates q^ are given by 
a solid curve. 



N=2 

N=4 

N =6 

-N,M 

0 

. 21 xl 0" 3 

.43xl0 “ 4 

. 20 x 10 

N,M 

E 1 

. 0207 

.0106 

.0055 

;n,m 

.18X10 -3 

.13xl0 -4 

.16X10 

N,M 

E 1 

.0213 

.0115 

.0060 

S N ' M 

.18xl0' 3 

.91xl0 -5 

.34x10 

e n,m 

.0048 

.0146 

.0010 


TABLE 4.1 
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Example 2 . We consider the same example as in Example 1 except that 

we used q*(x) = .15 + .10 tanh [20 (x- . 5) ] . Tests were conducted to 

compare the performance of our algorithms with quintic vs. cubic state 

approximations and cubic vs. linear parameter approximations. Typical 
''N M N.M 

values for $ and E ^ are given in Table 4.2 along with some 
typical graphs in Figures 4.2 and 4.3. In tabulating the results, 
we use the labels (N,M,S,P) with state approximations given by 

S = Q (quintic) or S - C (cubic) 
and parameter approximations given by 

P = C (cubic) or P = L (linear). 


(N,M,S,P) 

~N,M 

$ 

N,M 

-1— 

(8,7 ,Q,C) 

.546xl0“ 6 

.86xl0 -2 

(8,7,Q,L) 

.302xl0 -5 

.89xl0~ 2 

(8,7 ,C,C) 

.241xlo“ 6 

.12xl0 _1 

(8 ,7 ,C ,L) 

.224X10 -5 

.30xl0 -2 


TABLE 4.2 
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(12/7, C/C): N = 12/T1 =7/CUBIC STATES/CUBIC MWME1ERS 




(8,9, C/L): N=8/M=9/i 
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In Examples 3, 4 and 5 we consider a clamped-free or cantilevered 
beam with a variety of configurations at the free end, x=l. More 
precisely, in its most general form, we consider the beam with a 
rigidly attached tip body. The dynamics are given by the hybrid 
system of ordinary and partial differential equations (see [10] , [28] , 
[29], [30]) 

(4.5) pu tt + d2 ^ < 3^ d2u + q 2° 2u t^ = D ^ oDu ^ + 0<x< lr t>0, 

2 2 

(4.6) mu + mcDu - D{q D u + q D u } = - oDu + g, x=l, t>0, 

Li, L L JL ^ L 

2 2 2 

(4.7) mcu . + ( J + me ) Du + q. D u + q„D u. = - cODu + h, x = 1 , t>0, 

tt tt 1 2 t 

(4.8) u(t,0) = 0, Du ( t , 0 ) = 0, t>0, 

(4.9) u(0 ,x) = u Q (x) , u fc (0 ,x) = v Q ( x) , 0 < x < 1, 


where p=p(x) is the linear mass density of the beam, m is the mass of the 
tip body and J is its moment of inertia about its center of mass. The 
center of mass of the tip body is assumed to lie at a distance c from the 
tip of the beam directed along the tangent in the x-direction at the tip 
to the beam's neutral axis (see Figure 4.4 below). 

The second order term D{aDu} in (4.5) and the corresponding terms 
in (4.6) and (4.7) are the forces and moments which result from an axial 
loading a=o(t,x) . In the examples below, we only consider axial loads 



-3 
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which result from an acceleration a^ = a^ft) of base of the structure 
in the positive x-direction. In this case, we have (see [10], [29], [30]) 

a(t,x) = -a Q (t){m + f p(y)dy} . 

’'x 

The equations (4.6) and (4.7) represent respectively transverse and 
rotational equilibrium at the tip of the beam with g=g(t) denoting an 
externally applied load through the center of mass of the tip body in 
the transverse direction and h=h(t) an externally applied torque or 
moment . 

Note that setting p=l, m=J=c=0 and taking aQ=g=h=0 in equations 
(4.5), (4.6) and (4.7) leads to the dynamical equations for the standard 
cantilevered beam treated in Section 2. In the examples which follow we 
set 

f(t,x) = e sin2irt, 

took "true" values for the unknown parameters and used a quintic-spline 
based Ritz-Galerkin method to integrate the system (4.5) - (4.9) and 
generate displacement observations at times t^=.2i, i=0 , 1 ,2 , . . . , 5 at 
positions x^= . 2 5 j , j=2,3,4 along the span of the beam. The structure was 

assumed to be initially at rest (i.e., U q~ v q ==0 ^ • 

The modifications to the formulation of the approximation- schemes 
and the associated convergence theory necessitated by the presence of the 
tip appendages were briefly described and summarized in remarks in 
Sections 2 and 3 above. A complete and detailed discussion of our general 
approach in the context of inverse problems involving hybrid systems 
describing the vibration of beams with tip bodies can be found in [10] . 
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As might be apparent from the results of our numerical studies 
which will be presented below, clamped-free beams posed a somewhat 
stiffer challenge for our methods than did the simply-supported beams 
in the examples discussed previously. An inherent ill-posedness of 
the problem of estimating variable parameters in distributed systems 
was more evident here than when beams with simple boundary conditions 
were treated. 

An undesirable behavior (early onset of oscillations in the 
parameter estimates) may have resulted in part from the fact that in 
general in the presence of "higher order" boundary conditions it is 
more difficult to obtain accurate approximating solutions to the 
dynamical equations. 

Example 3 . In this example we simultaneously estimate a constant stiff- 

* 

ness coefficient, q^= .15 and a variable damping coefficient, 
q* (x) = . 01 (1. 5-tanh(20x - 10)) 

in a model of the form (4.5) - (4.9) for a cantilevered beam with no tip 
appendage (i.e. , with m = J = c = 0) . We took p (x) = 3 - x and = g= h = 0. 

Start up values for the optimization routine were chosen as q^ = .10 and 

q^Cx) = .015, 0 _< x <_ 1. 

Since in this example we are simultaneously estimating two parameters, 
the parameter space discretization index M is in fact a vector, M = 

(M^ , M^) of two indices with t-L corresponding to the discretization for 

q i , i = 1 , 2 . 

Using five cubic splines for the state approximation (N=4), two 
linear splines to discretize the sti-f fness coefficient (M^ = 1) and four 
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linear spline elements to discretize the damping coefficient (M^ = 3) 
we produced the estimates plotted in Figures 4.5 and 4.6. With eight 
cubic splines to discretize the damping coefficient all other 

approximation parameters left unchanged) we produced the estimates 
plotted in Figures 4.7 and 4.8. Using the labeling convention (N,M ^ , 
S,P ,P 2 ) we obtained 


(N,M L ,M ,S,P lf P 2 ) 

• ;n,m 

N,M 

E 1 

N,M 

E 2 

(4,1,3,C,L,L) 

.377 x 10 -6 

.908 x 10 -3 

.210 x lo' 2 

(4,l,5,C,L f C) 

-6 

.323 x 10 

.414 x 10 -2 

.256 x io -2 


The CPU times on the IBM 3081 for these two runs were 22.43 seconds and 


29.92 seconds respectively. 
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Example 4 . We consider a cantilevered beam with a point (c = J = 0) mass 

of magnitude m=1.5 rigidly attached to its free end. We took a^=l, 

■“ t -2t 

g(t) = 2e and h(t) =e . We simultaneously estimated a constant 

* 

stiffness coef fidient , q = .15 and a variable damping coefficient, 
q 2 (x) = .01(1.5 - tanh(3x- 1.5) ) . 

The start up values were taken to be = .1 and q 2 (x) = 0<_x<_l. 

with a five cubic spline based state approximation (N=4) , a two 
linear spline based discretization of the stiffness coefficient (M^ = 1) 
and a four linear spline (M 2 = 3) or a five cubic spline ( = 2 ) discre- 
tization of the damping coefficient we obtained the estimates plotted 
in Figures 4.9 - 4.12 with 


(n,m 1 ,M 2 ,s,p 1 ,p 2 ) 

cN,M 

?> 

N,M 

E 1 

N,M 

E 2 

(4,1,3,C,L,L) 

.747 x 10” 7 

.669 x 10 -3 

.364 x lo" 3 

(4,1,2,C,L,C) 

-7 

.404 x 10 

.572 x lo -3 

.702 x io“ 3 




0.0 0.5 1.0 

.'(4/L3, C/LL) : N=4/Fl 1 =l/I v l2 = 3/CUB I C STATES/LINEAR PARAMETERS/LINEAR PARAMETERS 

FIGURE A. 9 
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Example 5 . We simultaneously estimate a variable stiffness coefficient, 

q^(x) = . 1(.1.5 - tanh(20x - 10)) 

* 

and a constant damping coefficient = .015 for a cantilevered beam 
with tip body . We took m=1.5, c = . 1, J=.52, a = 1-, g(t) = 2e and 
h(t) =e 3t . start up values were chosen as q 3 (.x) = .15, 0<^x< 1, and q^ = 
. 01 . 

With a five cubic spline element state approximation (N=4) , a two 
linear spline based discretisation of the damping coefficient (M^ = 1) 
and either a four linear spline (,M^=3) or a nine cubic spline (M^ = 6) 
based discretization of the stiffness coefficient we obtained the results 
shown in Figures 4.13 - 4.16 with 


(n,m 1 ,m 2 ,s,p 1 ,p 2 ) 

;n,m 

N,M 

E 1 

N,M 

E 2 

(4,3,1,C,L,L) 

.695 x icf 7 

.174 x 10 -1 

.174 x 10 -3 

(4,6,1,C,C,L) 

.195 x io“ 7 

.181 x icf 1 

.550 x 10" 3 












0.5 1.0 

[ATES/CUBIC PARAMETERS/LINEAR PARAMETERS 

A. 15 
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5. Concluding Remarks 

We remark briefly on some of our findings as a result of the 

investigations described above as well as other related efforts. 

First, iti is well-known that "inverse" problems such as we have 

considered here are often ill-posed, lacking a certain "stability" 

(i.e., lacking a continuous dependence of the estimated parameters 
~N,M . 

q on the data {u. . j). This can often lead to serious difficulties 

13 

in computational efforts and it is sometimes helpful to use some type 
of regularization procedure in attempts to alleviate instabilities as 
well as speed up convergence . In much of our work we have taken an 
alternate approach, requiring that a compactness criterion be satisfied 
by the parameter set Q. Using arguments similar to some of those given 
in Section 2 above, one can argue that appropriate compactness assump- 
tions on the parameter sets will guarantee a type of stability (e.g. , 
see [1] ) . These compactness assumptions entail constraints that are 
sometimes "hard constraints" from a computational point of view, i.e. , 
it is desirable to implement them in order to obtain a well-behaved 
computational procedure. As we noted above, in some of our "test" 
calculations with beam models , we have needed to impose the pointwise 
constraints of Q (which are weaker than the compactness constraints 
in this particular problem) . As yet we have not tested the methods 
proposed in this paper with experimental data from beams (we are currently 
involved in projects to do this with beams of ordinary - e.g., steel- and 
composite materials) . However, we have had considerable experience 
using similar methods with experimental data in other problems (climatology , 
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bioturbation, insect dispersal) involving second order partial differential 
equation models. To date we have found the algorithms we propose well- 
behaved, often only requiring the imposition of parameter constraints that 

are easily implemented and incorporated into the optimization schemes. 

We believe that this will also be the case in dealing with data from 
experiments with beams and other elastic structures. 

Turning to a comment on the theoretical considerations of §2, we 
note that the presentations in [2] , [3] , [6] , [25] employ a semigroup approxi- 
mation approach (the Trotter-Kato theorem) that for the problems we 

2 

consider here would require H compactness of Q. The sesquilinear 

2 

form arguments of [19] appear to require H -weak compactness of Q as 
opposed to the ha compactness (in actuality, L compactness is sufficient 
and thus discontinuities can be allowed in the parameters if so desired) of 
this paper, [4] , [5] and [10] . Given our comments above regarding stability 
of the associated computational algorithms, this relaxation could well be 
of more than just theoretical interest. 

The methods described here did not, in general, perform as well when 
both a variable stiffness and damping coefficient were to be estimated 
simultaneously in our test example with beams. We are currently at work 
on a scheme that will deal more effectively with this more difficult 
class of problems. 

We note that both the theoretical and computational ideas outlined in 
this paper are applicable to a wide class of problems including beams 
modeled with the Timoshenko theories [3] and two dimensional elastic 
structures [8] . For further discussion of the advantages/disadvantages 
of weak vs. strong formulations in parameter estimation problems, we 


refer to the reader to [1] . 
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Finally, we remark that the numerical experiments outlined above 
(and others we have performed) are computationally intensive (even 
though for the examples of section 4 we always obtained convergence in 
40 or less iterations in the optimization algorithm).. Therefore, for 
conventional computational machines (e.g., sequential computers) use of 
these methods with experimental data could be expensive. However we 
believe that the approximation ideas and optimization schemes we are 
using offer great potential for use with emerging supercomputer technology 
(vector machines, attached array processors, parallel computers). In 
this regard, we are currently exploring the inherent parallelism and 
potential for vectorization of our algorithms and codes to develop fast, 
efficient software packages tailored to specific machine architectures. 

We are pursuing these efforts in the context of feedback control of 
distributed systems for elastic structures as well as for parameter 
estimation problems such as those considered above. Initial findings 
are quite encouraging and will be presented in a manuscript currently 
in preparation. 
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Appendix . 


We outline here the computation of the gradients via a costate 
formulation. We recall (.3.4) 


■ , N M m. 

$ (q) = j(u ;q 1 rq 2 ) 


V I N/*. M \ 

= l |u (t if x j? q ) - Uij 


1 ' D 


subject to (3.1) and (3.2) 


N N+l jj 

u (t,x) = y w w (t)B £ :(x) 
j=l ^ 3 


M 


M V .M. . - , 

q i = l q ij b j( X >' q = (q ij 5 
3=0 


N N NT 

where w = (w_ , . . , ,w_„ , _) satisfies (3,3): 

1 2N+2 


N*N N . M. N N 
Qw = K (q )w + f 


(A.l) 


N N . _ . N 
Q W (0) = W Q . 


N M M M 

We note that is independent of the parameters q = (q^»q 2 ) in 

the case considered here. We rewrite the cost criterion J as a 

functional over a continuum of t values by using Dirac functions 6 . 

For each N we have 


= l 


,T 




N+l 

I I 

k=l 


w^(t)B^(x.) ^u.(t)] 5 ( t-<-t . ) dt 
k k j 3 i 


(A. 2) 
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— M M M 2N+2 

where u.(t). = u... t. .. < t < t.. For q = (q, »q^!i and R 

3 13' 1-1 — 1 1 2 

vector functions w. and p, we define the associated Lagrangian 


m f T N+i 

(A. 3 ) ^Cw,p,q ). = y I y B, (x.)w (t). - u . ( t) ] 6 (t-t . ) dt 

i,j Jo kSl k 3 k 3 


T N» N M N 

P (t) [Q w(t). - K (q 2 w ( t) - f ]dt 


N 

Note that for any vector function p we have, for w = w the solution 
of (A. 1) , 


(A. 4) 


fj,. N, M, M 
$ (q) =iT(w (q ) ,p,q ) , 


so that 


(A. 5) 
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30"'" J<L ^r N M 3w" . , 3_zf. N M, 

= dizffw , p, q ] + r-f^(w ,p,q ) 


3q . . 
13 


3q. . 
13 


3q. . 

13 


for i = 1,2, j = 0,1,..., M, where is the differential of with 

respect to the vector function w. Thus the expression for the desired 

gradients can be greatly simplified if p is chosen so that djzftw N , p,q M 

for any variation function v (which we note must satisfy v(0) = 0 
N M 

since w Q is independent of q ). To. do this, we must compute the 
differential dSzf. We have 


dizf[w N ,p,q M ;v] = 


f T f 

{2 l ll 


B^(x .) w£(t)-u . (t) ] 5 (t-t.)]!’ b”(x.)v (t) 

o'- i,j k 3 3 i 3 


N, . . — 


,v] — 0 
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T N» N M 

+ p (t) [Q v(t) - K (q )v(.t)Udt 


T 

f fr- N T TN« NM 1 

= J |Iy Ct) 6(t-t i )v(t) +p i (t) [Q 1 v(t) - K (q M )v(t)]|dt 


N 2N+2 

where Y is the R vector function given by 


(A. 6) Y ^(t) 

X/ 


= < 


2 I rI B £ ( x )w£(.t) 
j k 


U (t)]B^(x ) , £ 
J J 


= 1,2, . . . ,N+1, 


0 £ = N+2, . , . , 2N+2 . 


Integrating by parts on the obvious, term in the expression for dJz? 7 
and using v(0) = 0, we find 


dLSf 
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lf*T N T NM r’NT 'l 

= J |-P Ct) Q - p (t)K (q“) + l Y M Ct) A 6(t-t i )jv(t)dt 


+ p T CT)Q N v(T) . 


Thus d Sf = 0 for all v if we choose p a solution of 


• T N T N M n N T 

p Ct) Q + p (t)K Cq ) - \ Y (t) 6(t-tJ = 0 
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P (T) = 0, 


.N . 


or since Q is a 2N+2 square symmetric matrix, 
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p(T) = 0. 
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Re turning to (A. 5) we thus find 


(A. 8) 
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N 

Recalling the definition of K Cg) - see (3.3) — and the representation 
(3.2), we find 


(A. 9) 
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for i = 1,2, j = 0,1,..,, M, where the 2N+2 square matrices 
are given by 


(A. 10) 
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, j — 0,1,. . . ,M 


with the N+l square matrix with elements 


(A. 12) 
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The costate equations (A, 7) can be transformed for convenience. Defining 


pet) = p<ti - iqV 1 


,t 
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0 i 1 
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Y" ( t i ) 
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we may compute the solution P of 


(A. 13) 


Q N P(.t) = -K N tq M ) T (p(t) + [Q N ] 1 l Y N (t )1 

l t.<t J 

N v N 

Q P CT) = 4 Y Ct.) . 
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This is then used to obtain p by 


(A. 14) pit) = Pit) + [Q N ] 1 l Y N Ct ) 

t.<t 


which in turn is used in (A. 9) to compute the desired gradients. 
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